Genome-wide identification of enhancers and transcription factors regulating the myogenic differentiation of bovine satellite cells

Background Satellite cells are the myogenic precursor cells in adult skeletal muscle. The objective of this study was to identify enhancers and transcription factors that regulate gene expression during the differentiation of bovine satellite cells into myotubes. Results Chromatin immunoprecipitation followed by deep sequencing (ChIP-seq) was performed to identify genomic regions where lysine 27 of H3 histone is acetylated (H3K27ac), i.e., active enhancers, from bovine satellite cells before and during differentiation into myotubes. A total of 19,027 and 47,669 H3K27ac-marked enhancers were consistently identified from two biological replicates of before- and during-differentiation bovine satellite cells, respectively. Of these enhancers, 5882 were specific to before-differentiation, 35,723 to during-differentiation, and 13,199 common to before- and during-differentiation bovine satellite cells. Whereas most of the before- or during-differentiation-specific H3K27ac-marked enhancers were located distally to the transcription start site, the enhancers common to before- and during-differentiation were located both distally and proximally to the transcription start site. The three sets of H3K27ac-marked enhancers were associated with functionally different genes and enriched with different transcription factor binding sites. Specifically, many of the H3K27ac-marked enhancers specific to during-differentiation bovine satellite cells were associated with genes involved in muscle structure and development, and were enriched with binding sites for the MyoD, AP-1, KLF, TEAD, and MEF2 families of transcription factors. A positive role was validated for Fos and FosB, two AP-1 family transcription factors, in the differentiation of bovine satellite cells into myotubes by siRNA-mediated knockdown. Conclusions Tens of thousands of H3K27ac-marked active enhancers have been identified from bovine satellite cells before or during differentiation. These enhancers contain binding sites not only for transcription factors whose role in satellite cell differentiation is well known but also for transcription factors whose role in satellite cell differentiation is unknown. These enhancers and transcription factors are valuable resources for understanding the complex mechanism that mediates gene expression during satellite cell differentiation. Because satellite cell differentiation is a key step in skeletal muscle growth, the enhancers, the transcription factors, and their target genes identified in this study are also valuable resources for identifying and interpreting skeletal muscle trait-associated DNA variants in cattle. Supplementary Information The online version contains supplementary material available at 10.1186/s12864-021-08224-7.


Background
Skeletal muscle is the largest tissue in the body and plays an important role in physiology [1,2]. Skeletal muscle from meat-producing animals is a major source of food for humans and animals. Adult skeletal muscle is composed of mostly muscle fibers [3]. A muscle fiber, also known as a myofiber, is a multinucleated muscle cell differentiated and fused from multiple mononuclear muscle cells called myoblasts. For most mammals, the total number of myofibers is determined prenatally [4]; thus, postnatal skeletal muscle growth results primarily from myofiber hypertrophy [4,5]. Postnatal myofiber hypertrophy, however, requires additional nuclei [6,7]. In postnatal animals, nuclei added to the existing myofibers are widely believed to come from satellite cells, which are mononuclear cells located near myofibers and are considered stem cells in adult skeletal muscle [8][9][10]. Satellite cells are normally quiescent but can be activated by muscle injury and nutritional and environmental changes [6,7]. Once activated, satellite cells become and proliferate as myoblasts and then differentiate and fuse with each other to generate new myotubes, the developing myofibers, or with existing myofibers to increase muscle fiber size.
A set of four transcription factors called myogenic regulatory factors (MRFs) are known to play important roles in myogenesis, the formation of muscle fibers from myoblasts or satellite cells [11,12]. These MRFs include myogenic differentiation 1 (MYOD1, also known as MyoD and MYF3), myogenic factor 5 (MYF5), myogenin (MYOG, also known as MYF4), and myogenic factor 6 (MYF6, also known as MRF4 and herculin). All four MRFs are specifically or preferentially expressed in skeletal muscle [12]. All four MRFs are basic helixloop-helix (bHLH) domain-containing transcription factors and regulate gene transcription by binding to the E-box sequence, CANNTG, where N is A, G, C, or T [13]. MYF5 and MYOD1 determine the myogenic lineage of stem cells in a redundant manner [11,12]. MYOG is essential for myoblast differentiation and fusion into myotubes [12,14]. MYF6 was thought to play a similar role to MYOG in myoblast differentiation, but a more recent study indicated an unexpected negative role of MYF6 in postnatal skeletal muscle growth [15].
Clearly, differentiation and fusion of myoblasts or satellite cells into myotubes is a key step in the development and growth of skeletal muscle. The objective of this study was to further understand the regulation of gene expression during the differentiation of bovine satellite cells into myotubes. Cattle are agriculturally important animals, and a better understanding of gene regulation during satellite cell differentiation could lead to the development of novel strategies to improve growth efficiency and meat quality in cattle. Enhancers are DNA sequences that enhance the transcription of associated genes when bound by sequence-specific transcription factors. Active enhancers are genomic regions that are bound by active transcription factors and that actively regulate gene transcription. Genomic regions containing active enhancers are found to be uniquely marked with H3K27ac, where lysine 27 of histone 3 protein is acetylated [16,17]. We began this study by identifying genomic regions with H3K27ac modification in bovine satellite cells before and during induced differentiation and fusion into myotubes through chromatin immunoprecipitation coupled with deep sequencing (ChIP-seq). To our knowledge, such an approach had not been taken to study the transcriptional mechanisms that control gene expression during the myogenic differentiation of bovine satellite cells or satellite cells from any species.

H3K27ac-marked enhancers in bovine satellite cells before and during differentiation
Four ChIP-seq libraries and two Input libraries constructed from bovine satellite cells immediately before and 2 days after induction of differentiation passed the quality control, and deep sequencing generated 23 to 40 million sequencing reads from these libraries (Table 1). Between 75 and 92% of these reads were uniquely mapped to the bovine genome, generating approximately 20 to 36 million uniquely mapped reads per library (Table 1).
Analyzing the uniquely mapped reads from each ChIPseq library against those from the corresponding Input library using the MACS peak calling program identified more than 30,000 and 50,000 H3K27ac-marked genomic regions, i.e., active enhancers, from before-differentiation (BD) and during-differentiation (DD) bovine satellite cells, respectively. A phantompeakqualtools analysis indicated that the four ChIP-seq libraries had normalized strand cross-correlation (NSC) values between 1.16 and 1.18 and relative strand cross-correlation (RSC) values between 0.98 and 1.06 (Additional File 1), indicating strong enrichment of reads in peaks [18].
A Pearson correlation analysis revealed that the H3K27ac-marked enhancer regions identified from two biological replicates, which corresponded to BD or DD satellite cells originally isolated from two different cattle, were highly correlated (Fig. 1A). A total of 19,027 H3K27ac-marked enhancers were consistently identified from two biological replicates of BD bovine satellite cells, while 47,669 H3K27ac-marked enhancers were consistently identified from two biological replicates of DD bovine satellite cells (Fig. 1B). A total of 5882 H3K27acmarked enhancers were found to be specific to BD bovine satellite cells, 35,723 H3K27ac-marked enhancers specific to DD bovine satellite cells, and 13,199 H3K27ac-marked enhancers common to both BD and DD bovine satellite cells (Fig. 1C, Additional Files 2, 3, 4).
Examples of H3K27ac-marked enhancers identified from BD and DD bovine satellite cells are shown in Table 2 and Fig. 2A. These enhancers were associated with the MYOG gene, which, as mentioned above, is a transcription factor essential for myoblast differentiation, and myosin heavy chain 3 (MYH3) gene, which, as indicated by its name, encodes a skeletal muscle-specific myosin heavy chain protein. As shown in Table 2   and Fig. 2A, MYOG-associated enhancers were marked with H3K27ac only in DD bovine satellite cells; MYH3associated enhancers were marked with H3K27ac in both BD and DD bovine satellite cells but more MYH3associated enhancers were marked with H3K27ac in DD than in BD bovine satellite cells. As shown by Fig. 2B, the increased H3K27ac modification to MYOG-and MYH3-associated genomic regions in DD satellite cells was accompanied by increased expression of both genes in these cells.

Genomic distribution of H3K27ac-marked enhancers in bovine satellite cells
H3K27ac-marked enhancers specific to BD or DD bovine satellite cells had different genomic distribution from the H3K27ac-marked enhancers common to both BD and DD bovine satellite cells. Whereas nearly 90% of H3K27ac-marked enhancers specific to BD or DD bovine satellite cells were located in the distal intergenic regions and introns, this percentage was only 60% for H3K27ac-marked enhancers common to BD and DD bovine satellite cells (Fig. 3A). Whereas approximately 6 and 17% of H3K27ac-marked enhancers specific to BD and DD bovine satellite cells, respectively, were located in the promoter regions, this percentage was almost 50% for H3K27ac-marked enhancers common to BD and DD bovine satellite cells (Fig. 3A). Whereas H3K27ac-marked enhancers specific to BD or DD bovine satellite cells were concentrated at 100,000 bp from the transcription start site (TSS), those common to BD and DD bovine satellite cells were concentrated at both 100 bp and 100,000 bp from the TSS ( Fig. 3B and C).

Expression levels of genes associated with and without H3K27ac modification in bovine satellite cells
We compared the expression levels of genes associated with H3K27ac modification with those without H3K27ac modification in bovine satellite cells. The transcriptome data of BD and DD bovine satellite cells were generated in a previous study [19]. As shown in Fig. 4A, in both BD and DD bovine satellite cells, genes associated with H3K27ac modification were expressed at levels nearly three-fold those of genes without H3K27ac modification. Overall, 5882 H3K27ac-marked enhancer regions Relative expression levels of MYOG and MHY3 mRNAs in BD and DD bovine satellite cells. *P < 0.05 (n = 4). Gene expression data was retrieved from a previous study [19] specific to BD bovine satellite cells were associated with 3064 protein-coding genes; 35,723 H3K27ac-marked enhancers specific to AD bovine satellite cells were associated with 7649 protein-coding genes; 13,199 H3K27acmarked enhancers common to BD and DD bovine satellite cells were associated with 6337 protein-coding genes. Obviously, many genes were associated with more than one H3K27ac-marked enhancers in bovine satellite cells. However, the numbers of H3K27ac-marked enhancers were not correlated with the expression levels of associated genes, regardless of the differentiation stage of the cells (Fig. 4B).

Functional terms enriched in genes associated with H3K27ac modification in bovine satellite cells
We performed gene ontology (GO) enrichment analyses on genes associated with H3K27ac modification in bovine satellite cells. Top biological processes and H3K27ac modification specifically in DD bovine satellite cells were related to skeletal muscle structure, development, and adaptation, and cell cycle arrest (Tables 3 and 4). It is interesting to note that many genes involved in melanosome and phagocytic vesicle were also associated with H3K27ac modification in DD bovine satellite cells ( Table 4). Most of the top molecular functions enriched in genes associated with H3K27ac modification specifically in DD bovine  satellite cells were related to growth factor binding and serine and threonine kinase signaling (Table 5).
Top biological processes, cellular components, and molecular functions enriched in genes associated with H3K27ac modification in BD bovine satellite cells included pentose metabolic process (Additional File 5), postsynaptic membrane (Additional File 6), and transmembrane receptor protein tyrosine kinase activity (Additional File 7), respectively, which were different from those enriched in genes associated with H3K27ac modification in DD satellite cells (Tables 3, 4 and 5).
Top biological processes, cellular components, and molecular functions enriched in genes associated with H3K27ac modification in both BD and DD bovine satellite cells included those related to proteasome and autophagosome (Additional Files 8, 9 and 10).

Transcription factor binding sites enriched in H3K27ac-marked enhancers in bovine satellite cells
Motif

Validation of the role of FOS and FOSB in bovine satellite cell differentiation
The motif enrichment analyses of H3K27ac-marked enhancers in BD and DD bovine satellite cells indicated that gene transcription during bovine satellite    Fig. 5A and B, in bovine satellite cells transfected with the negative control siRNA, the mRNA expression levels of these markers were significantly increased on day 3 of differentiation compared to their expression levels on the day before induction of differentiation. However, the mRNA expression levels of these markers in bovine satellite cells transfected with siRNA targeting FOS or FOSB mRNA were markedly lower than in cells transfected with the negative control siRNA. Morphologically, satellite cells transfected with siRNA targeting FOS or FOSB mRNA formed fewer and smaller myotubes than those transfected with the negative control siRNA (Fig. 5C). These data support a positive role for FOS and FOSB in bovine satellite cell differentiation and fusion into myotubes.

Discussion
Histone modification affects gene transcription by altering the accessibility of chromatin and the recruitment of transcription factors and cofactors to chromatin [20]. Large-scale mapping of histone modification has revealed that different types of enhancers are associated with different histone modifications. Active enhancers are associated with H3K27ac and histone 3 lysine 4 monomethylation (H3K4me1) modifications; primed or poised enhancers are marked with H3K4me1 but not H3K27ac; and silenced or repressed enhancers are often associated with H3K27me3 modification [21][22][23][24]. Based on these associations, ChIP-seq has been widely used to identify enhancers and other types of regulatory DNA regions in whole genomes [25][26][27][28][29][30][31]. In this study, we have identified 19,027 and 47,669 H3K27ac-marked enhancers in bovine satellite cells before and during differentiation, respectively. Identification of these enhances provides a valuable resource for understanding the mechanism that regulates gene expression during satellite cell differentiation, a key step in skeletal muscle development and growth. Compared to the consistent identification of 47,669 H3K27ac-marked active enhancers from two biological replicates of during-differentiation bovine satellite cells, only 19,027 H3K27ac-marked active enhancers were repeatedly identified from two samples of before-differentiation bovine satellite cells. This difference suggests much more active transcription factor binding to the genome, much more active recruitment of histone acetylases, and hence much more active H3K27ac modification in bovine satellite cells during differentiation than before differentiation. There is a possibility that this difference was caused by biological variation, as indicated by the large difference in the numbers of H3K27-marked enhancers identified from two samples of before-differentiation satellite cells. In our research on satellite cells, we have noticed satellite cells from different animals differ in differentiation potential in culture, and this difference suggests animal-to-animal variation in gene expression and histone modification in satellite cells.
Enhancers can be located upstream or downstream of TSS, in introns or exons, and near or distantly from the promoters [32][33][34]. Some enhancers can be located in intergenic regions several hundred kilobases away from TSS and control gene transcription by forming DNA loops with the promoters [35][36][37]. Genomic distribution analyses showed that H3K27ac-marked enhancers in bovine satellite cells have similar genomic distribution to enhancers in other types of cells or species [32][33][34]. However, the genomic location of H3K27ac-marked enhancers in bovine satellite cells varies with the differentiation stage of these cells. Whereas most of the differentiation stage-dependent H3K27ac-marked enhancers in bovine satellite cells were located in the distal intergenic regions, most of the differentiation stage-independent H3K27ac-marked enhancers in bovine satellite cells were located in the promoter regions. This study also showed that H3K27ac-marked enhancers specific to duringdifferentiation bovine satellite cells were associated with genes involved in muscle organization, adaptation, and development while H3K27ac-marked enhancers common to both before-and during-differentiation satellite cells were associated with genes involved in basic cellular functions and processes. These results suggest that distal enhancers are preferentially activated to increase the expression of genes determining the differentiation stage of satellite cells, whereas proximal enhancers or promoters are preferentially activated to increase the expression of genes maintaining the basic cellular function of satellite cells. This differentiation stage-dependent activation of distal and proximal enhancers in bovine satellite cells is apparently consistent with earlier findings that distal enhancers mediate the expression of cell type-or developmental stage-specific genes while core promoters and proximal enhancers are responsible for the expression of housekeeping genes [38,39].
This study showed that genes associated with H3K27ac modification were expressed at greater levels than those without H3K27ac modification in bovine satellite cells, regardless of the differentiation stage of the cells. This result supports H3K27ac as a histone marker for transcriptional activation [16]. This study also showed that many genes were associated with multiple H3K27ac-marked enhancers, but that the numbers of H3K27ac-marked enhancers were not correlated with the expression levels of associated genes in bovine satellite cells. These results suggest that multiple H3K27acmarked enhancers do not function in an additive manner to increase gene expression or that multiple H3K27acmarked enhancers are functionally redundant in bovine satellite cells. Indeed, recent studies using the CRISPR-Cas9 approach demonstrate that not every enhancer is functionally important and that most enhancers provide only a supportive or backup role in regulating gene expression [40][41][42].
H3K27ac modification at enhancers results from transcription factor binding and subsequent recruitment of histone acetyltransferases such as p300 and CBP [43,44]. In this study we have identified many transcription factors that may bind to H3K27ac-marked enhancers in bovine satellite cells before or during differentiation. Among the transcription factors that are predicted to bind to H3K27ac-marked enhancers in during-differentiation bovine satellite cells are MYOG and MYOD1, the MEF2 family transcription factors, the KLF family transcription factors, and the TEAD family transcription factors. Both MYOG and MYOD1 are known as the central transcriptional regulators of myoblast differentiation, and they regulate the expression of musclespecific genes by binding to the motif called E-box [45]. The MEF2 family transcription factors MEF2A, MEF2C and MEF2D [45,46], the KLF family transcription factors KLF3 and KLF5 [47,48], and the TEAD family transcription factors TEAD2 and TEAD4 [49][50][51] have also been shown to play a positive role in myoblast differentiation. Identification of the binding sites for MYOG, MYOD1, MEF2, KLF, and TEAD transcription factors among the top motifs enriched in H3K27ac-marked enhancers in during-differentiation bovine satellite cells validates the quality of active enhancers identified in this study.
This study shows that many other transcription factors regulate gene expression during satellite cell differentiation, and these other transcription factors include the AP-1 family of transcription factors (e.g., FOS). Enrichment of binding sites for the AP-1 family of transcription factors in active enhancers in during-differentiation bovine satellite cells is intriguing because the member of the AP-1 family of transcription factors JUN is known to antagonize the stimulatory effect of MYOD1 on myoblast differentiation [52]. However, overexpression of JUNB, a member of the AP-1 family of transcription factors closely related to JUN, increased hypertrophy and expression of the muscle-specific gene MYH4 in C2C12 myoblasts [53]. We have also validated a positive role of two members of the AP-1 family of transcription factors, namely FOS and FOSB, in driving bovine satellite cell differentiation in this study. Therefore, different members of the AP-1 family of transcription factors might have different effects on gene expression during myoblast or satellite cell differentiation.

Conclusions
In summary, we have identified tens of thousands of genomic regions associated with H3K27ac modification, i.e., active enhancers, in before-or during-differentiation bovine satellite cells. These enhancer regions contain binding sites for many transcription factors, including MYOG and MYOD1, whose role in myoblast or satellite cell differentiation is widely known, and the AP-1 transcription factors, AP-4, and many others, whose roles in myoblast or satellite cell differentiation are less known or unknown. These enhancers and transcription factors should be valuable for elucidating the mechanisms that mediate gene transcription during myoblast or satellite cell differentiation. Because myoblast or satellite cell differentiation is a key step of skeletal muscle development and growth, the enhancers, the transcription factors, and the genes targeted by these enhancers and transcription factors should be also valuable for identifying and interpreting skeletal muscle trait-associated DNA sequences and variants in cattle, which are agriculturally important animals.

Isolation and culture of bovine satellite cells
Skeletal muscle was collected from Angus-crossbred steers slaughtered at the Virginia Tech Meat Center. Satellite cells was isolated through pronase digestion and differential centrifugation as described before [54,55]. Satellite cells were cultured in growth medium for about a week before being induced to differentiate and fuse into myotubes. Differentiation of bovine satellite cells into myotubes was induced by replacing growth medium with differentiation medium. Growth medium consisted of Dulbecco's Modified Eagle Medium (DMEM), 10% fetal bovine serum (FBS) (R&D Systems, Minneapolis, MN, USA), 2 mM L-glutamine, and 1% Antibiotic-Antimycotic (100×) (ABAM). Differentiation medium consisted of DMEM, 2% horse serum (R&D Systems), 2 mM L-glutamine, and 1% ABAM. All cell culture was performed at 37 °C in a humidified, 5% CO 2 atmosphere. All cell culture reagents were purchased from ThermoFisher Scientific (Waltham, MA, USA) unless otherwise indicated.

ChIP assay
Satellite cells from two steers immediately before and 2 days after induction of differentiation were cross-linked in 1% formaldehyde for 10 min and then lysed in lysis buffer from the ChIP-IT kit (Active Motif, Carlsbad, CA, USA). Cell nuclei were suspended in ChIP buffer from the ChIP-IT kit and then sheared on ice by 10 pulses of 20-s sonication using a sonic dismembrator Model 100 at setting 3 (ThermoFisher Scientific) to generate chromatin fragments of 200 to 1000 bp. To identify the genomic regions associated with H3K27ac modification, chromatin fragments were incubated with an anti-histone H3K27ac antibody (ab4729, abcam, Cambridge, MA, USA) at 4 °C overnight with gentle rocking. The H3K27ac antibody-chromatin complexes were separated from unbound chromatin fragments using protein G-Dynal beads (ThermoFisher Scientific). Chromatin fragments immunoprecipitated by the H3K27ac antibody and those before immunoprecipitation (i.e., input chromatin) were reverse cross-linked by incubating them at 65 °C for 4 h. DNA was extracted and purified using spin columns from the ChIP-IT kit.

ChIP-seq library construction and sequencing
ChIP-seq libraries were prepared using the NEBNext ChIP-Seq Library Prep Reagent Set for Illumina (New England BioLabs, Ipswich, MA, USA), according to the supplier's instructions. Briefly, ChIP DNA was end repaired using T4 DNA polymerase, Klenow DNA polymerase, and T4 polynucleotide kinase. End-repaired DNA was then added with 3′ dA overhangs using exonuclease minus Klenow DNA polymerase and dATP. The dA-tailed DNA fragments were ligated to the sequencing adaptor. DNA fragments of approximately 300 bp were selected from the adaptor-ligated DNA using AMPure XP Beads. The size-selected DNA fragments were amplified in 12 cycles of PCR using an index primer and a universal PCR primer. Two ChIP-seq libraries were prepared from DNA immunoprecipitated from before-or duringdifferentiation bovine satellite cells originally isolated from two different cattle. Two Input-seq libraries were prepared from input DNA pooled equally from beforeand during-differentiation bovine satellite cells. Two biological replicates were used for ChIP-seq according to the ChIP-seq guidelines and practices proposed by the ENCODE and modENCODE consortia [18]. Each library was assessed for quality on a Bioanalyzer before being single-end sequenced on an Illumina Hiseq 2500 at the Genomics Sequencing Center at Virginia Tech.

Small interfering RNA (siRNA)-mediated gene knockdown
Bovine satellite cells in 12-well plates were transfected with 30 nM of siRNA targeting bovine FOS or FOSB mRNA using 6 μL of Lipofectamine RNAiMAX Reagent according to the supplier's instructions (ThermoFisher Scientific). A universal negative control siRNA (MIS-SION siRNA Universal Negative Control #1, Millipore Sigma, Burlington, MA, USA) was transfected as a negative control. The sense and antisense siRNA sequences targeting bovine FOS and FOSB mRNAs were CAG AAG AGA UGU CUG UGG CUU CUC U and AGA GAA GCC ACA GAC AUC UCU UCU G, and GAC AUG CCA GGA ACC AGU UAC UCC A and UGG AGU AAC UGG UUC CUG GCA UGU C, respectively. These siRNAs were confirmed to have at least 70% knockdown efficiency in pilot experiments. Following transfection, bovine satellite cells were cultured in differentiation medium for 48 h as descried above. The differentiation degree of satellite cells was assessed by quantifying mRNA expression of markers of differentiated myoblasts, including CKM, MYH2, MYH3, and MYOG, as previously described [19].

RNA extraction and reverse transcription-quantitative PCR (RT-qPCR)
Total RNA was extracted from bovine satellite cells using the Direct-zol RNA Miniprep Kit (Zymo Research, Irvine, CA, USA). Reverse transcription of total RNA into cDNA was performed using ImProm-II reverse transcriptase and random primers according to the manufacture's instruction (Promega, Madison, WI, USA). Quantitative PCR was performed using the SYBR Green chemistry as described previously [19]. The relative abundance of target mRNAs was calculated using the 2 -ΔΔCt method [67]. The Ct values for target mRNAs were normalized to the Ct values for HMBS, which was chosen as a reference gene because it was stably expressed in different conditions [68].

Statistical analysis
Gene expression data were analyzed by ANOVA. Two means were compared by t-test, and multiple means were compared by the Tukey test. All data are expressed as mean ± standard error.